*D-A03b  155  CONTROL  DATA  CORP  MINNEAPOLIS  MINN  DIGITAL  IMAGE  SYS ETC  F/G  8/2 

DIGITAL  CARTOGRAPHIC  STUDY  AND  BENCHMARK . (U) 

OCT  75  D J PANTON*  M E MUKPHY  DA AG53-75-C-0 195 

UNCLASSIFIED  ETL-0090  mi 


DIGITAL_^ARTOGRAPHIC^TUDY  AND ^ENCHMARK  m 
C FIRST  INTERIM  TECHNICAL  REPORT 


Prepared  for 

U.  8.  Army  Engineer  Topographic  Laboratories 
Fort-Rf lvoir.  Virginia 


DAAG53-75-C 


D.  J./ianton 
M.  E .j  Murphy 


Approved  For 
Distribution  Unlimited 


DIGITAL  IMAGE  SYSTEMS  DIVISION 
Control  Data  Corporation 
Minneapolis,  Minnesota 


processing  approach  and  a block  processing  approach.  Of  the  two  approaches. 
It  has  been  found  that  the  block  processing  approach  is  wore  adaptable  to  the 
requirements  of  digital  terrain  data  collection.  Therefore,  this  approach  Is 
being  Investigated  further  for  Implementation  In  an  array  of  fast,  micro- 
programmable  processors  to  provide  a benchmark  of  matching  system  parameters 
and  performance. 
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ABSTRACT 


The  key  problem  in  the  automatic  generation  of  digital  terrain 
data  is  the  matching  of  conjugate  points  on  a stereo  pair  of  aerial 

photographs  in  an  accurate  and  timely  manner.  The-ewh  deoowibed 
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4«i  this  report  fcentera  around /the  development  of  suitable  algorithms 
and  systems  procedures  to  perform  this  matching  task  In  an  all-digital 
environment . 
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Two  approaches  to  automatic  image  matching  -have  baifo/lnvcsti- 
gated;--a  strip  processing  approach  and  a block  processing  approach. 
Of  the  two  approaches,  it  has  beer.> found  that  the  block  processing 

approach  is  more  adaptable  to  the  requirements  of  digital  terrain 
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data  collection^  Therefore,  -this  approach A is  being  investigated 
further  for  implementation  in  an  array  of  fast,  microprogrammable 
processors  to  provide  a benchmark  of  matching  system  parameters  and 
performance . 
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1.0  INTRODUCTION 


Automatic  systems  for  the  compilation  of  digital  terrain  data  must  be  fast 
accurate  and  flexible.  In  addition,  the  cost  of  such  systems  must  be  in  pro- 
portion to  the  quality  of  data  they  generate.  The  need  for  automatic  carto- 
graphic systems  is  unquestioned,  considering  the  volume  of  photographic  data  in 
need  of  processing.  What  is  Important  to  consider  is  the  optimum  technology 
required  to  perform  this  task  of  cartographic  data  processing. 


This  technical  report  is  a summary  of  the  .first  phase  of  ongoing  work 
in  the  development  of  fast,  accurate  and  flexible  algorithms  and  systems  to 
produce  the  desired  mapping  products. 


The  primary  area  of  consideration  in  generating  digital  terrain  data  is 
the  automatic  matching  of  conjugate  points  in  a stereo  pair  of  aerial  photographs 
This  is  the  area  in  which  computational  speed  is  of  utmost  importance  and  in 
which  a great  deal  of  algorithm  sophistication  is  required.  Clearly,  all  of  the 
research  is  not  yet  in  regarding  the  optimum  approach  for  such  a task. 


Therefore,  two  approaches  to  the  problem  of  image  matching  which  have 
been  used  in  the  past  by  Control  Data  for  automatic  registration  and  change 
detection  purposes  are  being  investigated.  The  study  described  in  this 
report  involves  the  redesign  and  modification  of  the  two  approaches  to 
efficiently  and  accurately  process  stereo  imagery.  The  overall  goal  of  the 
effort  is  to  choose  the  better  approach  for  Implementation  in  fast,  micro- 
programmable  processors  as  a demonstration  and  verification  of  feasibility. 


The  two  image  matching  approaches,  called  the  strip  approach  and  block 
approach,  have  been  generally  described  in  a previous  technical  report  [l]  , 
and  these  aspects  will  not  be  re-discussed  here.  Only  the  modifications,  new 
aspects,  and  new  results  will  be  treated.  Suffice  it  to  say  that  since 
both  approaches  utilize  the  correlation  coefficient  as  the  image  matching 
metric,  what  the  study  is  really  comparing  are  the  two  philosophies  behind  the 
matching  approaches. 


The  strip  processing  approach  is  very  fast  and  production  oriented.  It  is 
fast  because  it  incorporates  a line-by-line  error  correcting  process.  That  is, 
the  process  knows  where  on  the  image  to  move  next  based  on  the  errors  it  sees 
now.  On  the  other  hand,  the  block  approach  is  a slower  process.  Incorporating  a 
a rather  deterministic  predict-ahead  mechanism.  The  study  has  shown  that 
because  of  these  differences  in  approach,  differences  in  the  output  digital 
terrain  data  occur  and  a number  of  trade-offs  emerge. 


2.0  MODIFICATION  OF  STRIP  APPROACH:  INTRODUCTION 

Throughout  this  section  of  the  report,  we  will  be  addressing  the  problem 
of  successfully  matching  stereo  imagery  with  an  automatic  strip  process.  This 
automatic  strip  processor  will  be  referred  to  as  program  TRAK.  We  will  discuss 
enhancements  to  TRAK  made  possible  by  the  nature  of  stereo  imagery.  In 
addition,  we  will  also  address  problems  identified  by  an  earlier  attempt  to 
use  the  strip  process  for  this  application.  We  will  discuss  solutions 
implemented  for  those  problems  and  will  indicate  results  both  visual  and 
statistical  which  we  were  able  to  achieve. 

2.1  INCORPORATION  OF  EPIPOIAR  Y-AXIS  CONTROL 

Previous  applications  of  the  strip  process,  in  particular,  side  looking 
radar,  ERTS  multispectral  data,  and  high  resolution  aerial  photography,  have 
Included  a correlation  process  using  5 patch  sites  separated  by  a distance 
of  from  1 to  7 cells  as  indicated  in  Figure  2-1. 

The  effective  size  of  the  5 sites  is  some  patch  height  (user  defined) 
by  some  patch  width.  A weighted  accumulation  of  the  required  vector  dot  pro- 
ducts provides  an  effective  patch  width  dimension  thereby  reducing  required 
memory  for  the  actual  inner  products  necessary  -to  compute  correlation.  Using 
these  5 correlations,  an  interpolated  maximum  point  on  the  correlation  surface 
is  found  as  illustrated  in  Figure  2-2.  Deviations,  or  spatial  error  terms, 
are  computed  and  then  incorporated  into  warp  equations  as  shown  in  Figure  2-3. 
Subsequently,  these  error  terms  drive  the  correction  process. 

• For  the  application  of  matching  stereo  imagery,  north  and  south  sites 
were  eliminated  from  the  correlation  and  subsequent  warp  update  process 
Epipolar  line  geometry  was  incorporated  into  the  process  thereby  providing 
y-axis  control.  X-axis  matching  is  as  performed  in  previous  applications. 
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Figute  2-3.  Relationship  of  the  Maximum  Correlation  Coefficient 
Location  as  Used  to  Update  Prediction  Equation. 
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2.2  PROBLEMS  IDENTIFIED  BY  PREVIOUS  EXPERIENCE 


Based  on  a previous  study  done  by  Digital  Image  Systems  Division  for  USAETL, 
the  following  problems  were  identified  in  the  strip  process:  Fi. jt,  the  prob- 

lem of  effectively  modeling  the  terrain  in  the  creation  of  a synthetic  scan 
line  for  subsequent  correlation  and  match  update;  second,  the  problem  of 
allowing  the  warp  update  process  enough  freedom  to  follow  terrain,  and  yet 
retain  enough  stability  in  areas  of  dissimilarity;  and  third,  the  problem 
of  interlocking  strip  information  to  effectively  guide  strips  in  hard  to 
match  areas.  In  the  following  3 sections,  we  shall  address  each  of  these 
problems  and  the  solutions  we  implemented  for  each. 

2.2.1  APPROXIMATING  HIGH  ORDER  EQUATION  IN  OBTAINING  SYNTHETIC  SCAN  LINE 

It  is  well  known  that  the  successful  matching  of  stereo  images  requires 
the  removal  of  x-parallax  or  differences  in  the  placement  of  matching  pixels 
in  stereo  images  due  to  terrain  relief.  Since  the  TRAK  program  is  a line- 
by-line  process,  it  requires  the  creation  of  a synthetic  scan  line  from  the 
dependent  image  of  the  stereo  pair  for  each  line  of  the  independent  image. 

The  proper  creation  of  this  synthetic  scan  line  is  of  utmost  importance  in 
determining  the  current  parallax  signal,  which  in  turn  is  used  to  drive  the 
process.  The  building  of  this  equivalent  synthetic  line  is  accomplished  by 
defining  positions  within  a buffer  containing  a section  of  the  dependent 
image  which  are  the  equivalent  of  the  strip  centers  in  the  independent  image. 
This  relationship  is  shown  in  Figure  2-4.  Interpolation  between  strip 
centers  is  linear  and  uses  a nearest  neighbor  criterion  for  the  actual  access  of 
matching  pixels.  By  retaining  the  linear  interpolation  technique  and  moving 
Atrip  centers  closer  together,  a much  better  approximation  to  the  matching 
synthetic  scan  lino,  which  is  really  a function  of  the  terrain  being  imaged, 
can  be  obtained.  This  effectively  gives  TRAK  the  capacity  to  closely  approxi- 
mate a warp  equation  of  an  order  which  is  entirely  dependent  on  the  terrain 
being  encountered.  In  addition,  correlation  patch  shaping  is  effected  each 
line  since  inner  product  memories  can  be  updated  with  information  resulting 
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from  the  previous  line's  parallax  signal. 
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2.2.2  CONSTRAINT  VERSUS  FLEXIBILITY  IN  THE  CHOICE  OF  SYSTEM  DAMPING 

System  damping  in  the  TRAK  program  may  be  defined  as  the  amount  of 
filtering  or  smoothing  needed  to  sufficiently  remove  noise  effects  in  the 
error  signal  values  produced  by  the  correlation  process.  There  are  two  unique 
and  yet  related  methods  by  which  the  present  version  of  TRAK  is  damped;  the 
first  is  the  rate  at  which  inner  products  used  for  the  correlation  process 
are  allowed  to  accumulate  (definition  of  effective  patch  width) , and  the 
second  is  the  rate  at  which  error  signals,  resulting  from  location  of  the 
interpolated  maximum  on  the  correlation  surface,  are  allowed  to  influence  the 
predicted  warp  positions.  In  previous  applications,  for  which  the  strip  process 
has  been  used,  studies  showed  a patch  width  of  10  lines,  providing  a correla- 
tion coefficient  based  on  a combination  of  data  smoothed  over  approximately 
7.5  lines  and  providing  independent  spatial  error  signals  approximately 
every  20  lines,  was  sufficient  for  most  applications.  It  was  only  necessary 
to  control  the  rate  at  which  these  error  signals  were  allowed  to  update 
matching  positions  in  the  dependent  image  to  achieve  satisfactory  results. 

For  this  application, however,  it  was  discovered  the  frequency  of  independence 
in  the  error  signals  was  required  at  a rate  more  often  than  20  lines.  In 
addition,  it  was  found  a much  larger  percentage  of  the  error  signal  had  to 
be  incorporated  during  the  warp  update  computations  to  maintain  proper  regis- 
tration in  areas  of  rapidly  changing  terrain.  In  obtaining  the  results  in- 
dicated later  in  this  report,  an  effective  patch  size  of  length  20  pixels  and 
width  defined  by  a smoothing  factor  with  half-decay  time  of  5 lines,  coupled 
with  an  error  signal  smoothing  factor  of  6 lines  were  used.  This  combination 
appears  to  have  provided  needed  flexibility  in  adjusting  position  changes  due  to 
terrain  encountered. 
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2.2.3  METHOD  OF  INTERLOCKING  STRIP  INFORMATION  TO  GUIDE  STRIPS  IN  HARD  TO 

MATCH  AREAS 

Throughout  this  discussion,  we  will  refer  to  the  task  of  guiding  strips 
In  hard  to  match  areas  as  harnessing.  In  previous  applications  of  the  strip 
process,  the  technique  used  in  satisfying  this  requirements  is  as  Illustrated 
in  Figure  2.5.  As  the  figure  suggests,  current  strip  positions  are  weighted 
by  some  measure  of  their  accuracy  and  a first  to  fourth  order  equation  (user 
defined)  is  derived  using  a least  squares  technique.  Current  strip  positions 
are  compared  with  least  squares  estimates  and  any  strips  exceeding  some  cir- 
cular tolerance  limit  are  corrected  by  averaging  the  current  position  and 
least  squares  estimate.  Since  the  harness  equation  can  be  first  to  fourth 
order  in  complexity,  this  technique  has  proven  acceptable  in  modeling  all 
types  of  distortion  encountered  to  date  in  applications  for  which  the  strip 
process  has  been  used. 

In  the  case  of  stereo  images,  this  technique  proved  unacceptable.  The 
problem  becomes  one  of  successfully  approximating  a frontal  curve  defined  by 
current  strip  positions.  The  order  of  the  equation  necessary  for  this  task 
is  a function  of  the  terrain  which  is  currently  being  matched.  It  becomes 
apparent  as  the  length  of  the  scan  line  being  matched  increases,  and  the  terrain 
within  that  scan  lines  varies,  the  complexity  of  the  equation  necessary  to 
model  that  scan  line  becomes  increasingly  great. 

Essentially  what  we  implemented  in  order  to  overcome  this  problem  was  a 
number  of  piecewise  linear  approximations  to  the  current  frontal  curve  as 
illustrated  in  Figure  2.6.  It  was  found  that  the  information  from  6 strips 
provided  the  needed  stability  for  the  least  squares  estimate  to  control  strips 
within  any  particular  segment  of  the  frontal  curve  within  some  circular  toler- 
ance, typically  3 to  5 pixels.  It  was  our  conclusion  that  this  would  rarely 
be  the  case.  In  addition,  the  piecewise  linear  approximation  has  the  added 
advantage  of  allowing  much  greater  scan  line  lengths  to  be  processed  without 
the  frontal  curve  modeling  problem  Inherent  in  the  greater  length. 
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Figure  2-5.  HARNESS  TECHNIQUE  FOR  SIDE  LOOKING  RADAR 
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FIGURE  2-6.  HARNESS  TECHNIQUE  USED  FOR  MATCHING  STEREO  IMAGES 
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3.0  MODIFICATION  OF  BLOCK  APPROACH 


Under  • previous  contract,  a block  processing  software  system  was  de- 
signed and  delivered  to  USAETL  for  use  with  the  DIMES  System  [l] . This 
block  system  was  research"or tented  and  capable  of  being  modified  in  a 
variety  of  ways  for  different  applications. 

For  the  present  image  matching  effort,  considerable  changes  have  been 
made  to  this  system,  but  the  underlying  philosophy  of  block  processing  or 
area  correlation  has  remained  constant.  Under  this  concept  a correlation 
area  of  certain  dimensions  is  defined  on  one  image  of  the  stereo  pair  and  a 
correlation  search  for  the  matching  area  on  the  other  image  is  initiated. 
This  section  describes  the  particular  algorithms  that  have  been  substituted 
or  added  under  the  current  work  effort. 

3.1  DATA  MANAGEMENT 

Previously  the  block- processor  was  equipped  with  a data  manager  that 
would  allow  correlation  to  take  place  at  random  anywhere  in  the  imagery. 

The  cost  of  this  flexibility  was  a reduction  in  speed. 

Therefore,  in  the  present  effort  the  assumption  was  made  that  the 
imagery  is  to  be  traversed  in  an  orderly  manner  from  one  end  to  the  other, 
that  is,  from  scan  line  1 to  scan  line  N.  Data  management,  then,  is  greatly 
simplified,  consisting  of  the  maintenance  of  two  buffer  windows,  one  on  each 
image  of  the  stereo  pair,  that  move  across  the  image  with  the  correlation 
process . 

i 

To  illustrate  this  data  manegement,  reference  is  made  to  Figure  3-1. 

The  buffer  window  on  the  independent  image,  hereafter  to  be  called  Image  A, 
la  wide  enough  to  include  the  number  of  scan  lines  required  for  a predeter- 
mined patch  width.  The  buffer  window  on  the  dependent  image,  to  be  called 
Image  B,  contains  as  many  scan  lines  as  is  required  for  the  patch  size 
search  area,  and  the  estimated  warp  between  Image  B and  Image  A. 
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Hew  scan  lines  are  read  into  the  buffers  only  when  needed  as  the  corre 
la t ion  patches  progress  across  the  image.  Image  matching  occurs  first  for 
all  patches  In  Image  A whose  centers  lie  on  the  same  scan  lines  before  the 
patches  move  on  to  subsequent  scan  lines  and  before  new  data  is  added  to 
the  buffers. 


Direction 
of  Processing 


Scan 

Line 

Direction 


•Correlation 
Search  Area 


corre  latio; 
patch 


Image  A 
Buffer  Window 


Image  B 
Buffer  Window 


Figure  3-1  Data  Management  Scheme 


Regarding  the  Image  B buffer  window,  there  are  three  distinct  cases 
of  buffer  management  that  can  occur  for  every  match  point  determination. 
Figure  3-2  depicts  these  cases. 


Ihe  optimum  situation  exists  when  the  buffer  window  is  wide  enough 
such  that  case  1 does  not  occur. 


Case  1 - Patch  and  search  area  behind 

current  buffer  window  placement* 
either  enlarge  buffer  or 
back  up  the  data 


Case  2 - Patch  and  search  area 
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no  buffer  management  necessary 


-*»  Case  3 - Patch  and  search  area  ahead  of 
{ current  buffer  window  placement 

i move  window  ahead  by  reading 

“ in  more  scan  lines 


Figure  3-2  Three  Cases  of  Image  B Buffer  Management 
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3.2  CORRELATION  ALGORITHM 


In  the  determination  of  a match  point,  it  is  necessary  to  place  the 
correlation  patch  center  at  each  site  of  the  Image  B search  area  and 
accumulate  over  the  patch  area  the  necessary  gray  scale  sums  and  cross 
products  needed  for  the  computation  of  the  correlation  coefficient.  The 
net  result  is  a value  of  the  correlation  coefficient  for  each  site  of  the 
search  area. 

Previously  in  block  matching  systems,  the  computation  of  sums,  cross 
products,  and  correlation  coefficients  has  occurred  in  total  fer  each 
distinct  patch  placement  in  the  search  area.  In  this  scheme  there  are  a 
great  deal  of  redundant  computations,  and  each  Image  B pixel  is  accessed 
repeatedly,  once  for  each  patch  placement  that  it  is  contained  in. 

Under  the  current  matching  effort,  the  basic  philosophy  behind  the 
correlation  algorithm  is  to  access  each  pixel  only  once  for  a given  search 
and  to  accumulate  its  gray-scale  value  when  it  is  available  in  all  the  sums 
and  cross  products  for  which  it  has  influence.  The  actual  correlation 
algorithm  directs  patch  placement  over  a search  area  that  extends  in  only 
one  dimension  around  a predicted  point  on  Image  B.  Thus,  correlation  values 
are  generated  on  a line  segment  coincident  with  the  eplpolar  line  passing 
through  the  predicted  point. 

The  mechanics  of  the  algorithm  are  illustrated  in  Figure  3-3.  As  each 
pixel  of  the  A image  patch  is  accessed,  its  value  is  accumulated  in  the  patch 
eta  and  sum  of  squares.  Likewise,  as  each  pixel  of  the  B image  search  area 
la  accessed,  its  value  is  accumulated  in  its  respective  column  sum  and  column 
sum  of  squares.  Also,  cross  products  are  generated  for  each  patch  placement 
along  the  Image  B search  segment.  When  the  accumulation  of  sums  and  cross 
products  is  complete,  the  variance  of  the  Image  A patch  is  computed  and 
the  array  of  column  sums  and  the  array  of  column  sums  of  squares  are 
traversed  to  generate  an  Image  B variance  and  a covariance  for  each  site  of 
the  search  segment.  The  array  traversal  involves  the  addition  of  the  next 
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column  sum  and  the  subtraction  of  the  last  column  sum  from  a running  total 
to  simulate  the  movement  of  the  patch  along  the  search  segment. 


Image  A 
Patch 


J 

Patch  Sum 

Patch  Sum 
of  Squares 


Image  B 
Search  Area 


Figure  3-3  Correlation  Algorithm  Detail 


The  final  result  of  the  algorithm  is  a.,  array  of  correlation  coefficients 
of  the  form 

• . ,m  (A.  & 

>J  var(A)  var(B)  , 

one  value  for  each  placement  of  the  patch  along  the  search  segment. 

3.2.1  CORRELATION  MAXIMUM  DETERMINATION 

It  has  been  shown  in  previous  studies  that  it  is  not  sufficiently 
accurate  for  cartographic  purposes  to  merely  find  the  maximum  value  of  the 
correlation  coefficient  along  the  search  segment  at  a pixel  center. 
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Rather  it  is  necessary  to  interpolate  the  correlation  maximum  to  a fraction 
of  a pixel.  , 

Referring  to  Figure  3-4,  the  correlation  maximum  at  a pixel  center  is 
found  and  then  by  using  one  value  on  either  side  of  this  maximum,  a parabolic 
curve  Is  fit  over  the  three  values.  The  point  along  the  search  line  segment 
at  tdiich  the  derivative  of  the  parabolic  curve  is  zero,  is  the  point  of  maximum 
correlation. 


Figure  3-4  Correlation  Maximum 

A problem  exists  when  the  maximum  correlation  value  at  a pixel  center 
occurs  at -either  extremity  of  the  search  segment.  In  this  case,  a parabolic 
fit  cannot  be  performed,  and  the  algorithm  is  designed  to  flag  this  situation 
and  return  the  pixel  center  as  the  point  of  maximum  correlation.  The  under- 
lying idea  here  is  that  this  match  point  is  suspect  and  a candidate  for 
further  processing. 

3.2.2  MATCH  POINT  DETERMINATION 

The  result  of  one  loop  through  the  correlation  algorithm  is  a pair  of 
conjugate  match  points  of  the  form  (x,y,u,v)  where  x,y  are  the  coordinates 
of  the  point  on  Image  A and  u,v  are  the  coordinates  on  Image  B.  These  points 
■ay  be  expressed  by  the  algorithm  in  either  digital  scan  coordinates  or 


photo  coordinates 


Since  the  correlation  algorithm  performs  the  correlation  search  in 
only  one  direction,  that  is,  along  an  epipolar  line  segement,  what  is  really 
searched  and  found  by  gray-scale  correlation  is  the  u coordinate.  The  v 
coordinate  of  the  match  point  is  computed  directly  from  epipolar  geometry 
parameters;  thus  eliminating  the  need  for  a correlation  search  in  the  v 
direction. 

3.3  PREDICTION  SCHEME 

Prediction  in  automatic  matching  systems  involves  the  accurate  acquisi- 
tion of  the  next  point  to  be  used  to  center  a search  area  based  on  previously 
matched  points.  In  the  previous  block  matching  systems  such  prediction  was 
made  using  a large  number  of  match  points  (in  the  neighborhood  of  20  to  30) . 
Minimal  use  of  stereo  geometry  was  used  in  the  prediction,  snd  because  of 
the  large  number  of  points  used,  the  prediction  was  rather  global  in  nature. 
This  required  the  use  of  large  correlation  search  areas  to  compensate  for 
the  local  image  deviations  from  the  global  prediction.  However,  it  was 
found  that  large,  two-dimensional  search  areas  contributed  to  the  instability 
of  the  prediction  mechanism  in  hard-to-correlate  areas. 

Iherefore,  it  was  advisable  in  the  current  matching  effort  to  design  a 
prediction  scheme  that: 

• uses  as  much  stereo  geometry  as  possible 

• is  more  locally  valid 

• is  accurate  enough  to  reduce  the  search  area  to  a minimum 

• relies  less  heavily  on  the  correlation  coefficient 

Such  a predictor  was  implemented  and  is  described  in  the  following 
sections . 
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3.3.1  EPIPOLAR  CONTROL 


Ab  was  mentioned  above,  the  correlation  search  for  the  v coordinate 
in  the  former  matching  system  was  replaced  by  the  direct  computation  of  the 
v coordinate  from  epipolar  geometry.  Referring  to  Figure  3-5,  the  large 
dots  represent  previous  matched,  conjugate  points. 

Image  A Image  B 


Figure  3-5  Prediction  Mechanism 


The  x's  within  the  patches  indicate  the  points  to  be  matched  next.  The 
sequence  of  matching  over  the  Image  A regular  grid  proceeds  from  point  i,j 
to  i,j+l  to  i,j+2,  etc.  until  one  line  of  blocks  is  complete.  Then  matching 
continues  with  line  i+1. 

Using  the  photo  coordinates  of  the  center  point  of  the  patch  i+1,  j+1 
on  the  A image  and  the  relative  orientation  elements  of  Image  B to  Image  A, 
the  equations  of  corresponding  epipolar  lines  are  determined.  Knowledge  of 
where  the  epipolar  line  lies  on  the  B Image  essentially  determines  the  v 
coordinate  of  the  match  point  without  need  for  a correlation  search  in  the 
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v direction.  What  is  necessary  is  an  accurate  prediction  of  the  u coordin- 
ate along  the  epipolar  line. 


3.3.2  RATE  OF  CHANGE  FUNCTIONS 


The  scheme  used  for  predicting  the  u coordinate  is  based  on  tp! £pt, 
the  rate  of  change  of  the  Image  B u coordinate  with  respect  to  the  Image  A 
x coordinate.  This  function  is  kept  and  updated  independently  for  each 
path  of  blocks,  j,  running  along  the  images  in  the  parallax  direction. 

/Sx  is  the  predetermined,  regular  interval  between  match  points  on  Image  A and 
Ai's  are  computed  for  each  path  j as  new  match  points  are  found.  In 
Figure  3-5,  for  example. 


To  obtain  a prediction  for  the  u coordinate  for  match  point  pair 
(i+l,  j+1) , a combination  of  the  neighboring  values  of  ip/tx  may  be  used. 
For  the  present  image  matching  effort,  the  following  prediction  formula 
was  found  most  successful  for  the  mountainous  areas  of  the  Phoenix  imagery 


where  ui+1  is  a prediction  to  be  refine*1  by  the  correlation  search. 
When  the  correlation  algorithm  has  found  the  tree  match  point  u^+^  j+^ 
along  the  epipolar  line,  then  **  aP**ted  using  the  new 

match  point  for  that  path  and  the  predict  ion  Mchanism  moves  on  to  match 
location  (i+l,j+2). 


The  advantage  of  such  a prediction  scheme  is  that  it  can  be  as  local 
or  as  global  as  desired,  depending  upon  how  many  and  the  placement  of 


neighboring  rate  of  change  functions  that  are  used  In  a prediction.  Moreover 
the  fij/flc  function  is  proportional  to  the  actual  terrain  slope  as  it 
appears  in  terms  of  Image  A.  The  following  formulation  and  its  derivations 
show  this  relationship: 


where  is  the  terrain  slope  with  respect  to  the  x coordinate,  B is  the 
baseline  distance  between  exposure  stations,  f is  the  focal  length  of  the 
cameras  and  P^  and  are  parallax  values  over  the  interval  in  which 

the  changes  are  computed.  Tills  equation  holds  only  for  truly  vertical  photo 
photographs  and  is  set  down  here  to  point  out  that  £u//sc  varies  inversely 
with  the  terrain  slope  and  is  close  to  1.0  in  flat  terrain.  As  depicted 
in  Figure  3-6,  &i/&c  takes  on  values  different  from  1.0  depending  upon 
whether  the  terrain  is  ascending  or  descending  in  the  process  direction. 


Figure  3-6  Effect  of  Terrain  Slope  on  flu/ /sc  Function 


It  has  been  found  that  for  stereo  base/height  ratios  in  the 


neighborhood  of  0.6  and  for  terrain  slopes  that  do  not  exceed  ± 45  degrees 


the  rate  of  change  function  is  generally  within  the  range 


In  terms  of  automatic  matching,  the  stability  of  a prediction 
mechanism  using  the  £m//?c  function  is  very  much  dependent  on  both  the 
choice  of  /sc  and  the  size  of  the  B Image  search  segment  in  scan  lines 


The  rationale  behind  the  optimum  choice  of  these  parameters  is  to  obtain  a 
prediction  within  the  allowable  limits  of  the  W & function  so  that 
correlation  is  performed  in  but  a minimum  area.  In  tibls  way,  less 
emphasis  is  placed  on  the  correlation  coefficient  finding  a reliable 
match  point  over  an  extended  area.  Table  3-1  hypothesizes  a 7-line  search 
segment  centered  at  a predicted  point  and  summarizes  the  behavior  of  the 
fli/ /pc  function  for  given  values  of  1st  and  for  a correlation  maximum 
occurring  at  each  site  of  the  search  segment. 


Table  3-1  THE  EFFECT  OF  SEARCH  SEGMENT  SIZE  ON  THE  tn/ &■ 
PREDICTING  FUNCTION  ' 


i i i ■ • ■ 

OOOXOOO- 

-3  -2  -1  0 +1  +2  +3: 


tti  (In  Lines  from  Predicted  Point) 


1.1083 

1.166 

1.25 

l.<062 

1.125 

1.187 

The  table  is  set  up  for  the  case  where  the  center  of  the  search  segment  is 
.at  £u*flx,  the  case  of  flat  terrain  starting  an  inflection  upward  or  down-  ' 
ward . 

The  dashed  area  of  the  table  indicates  the  search  segment  lengths  that 


are  critical  for  maximum  stability  of  prediction  for  various  values  of  &c. 
That  is,  for  to c values  of  2,  4 and  8,  the  optimum  search  segment  lengths 
are  3,  5 and  7 respectively.  What  this  means,  for  example,  is  that  if  a 
7- line  search  segment  were  used  with  a Ax  of  2 lines  and  further,  if  the 
area  being  searched  is  lacking  in  feature  or  grossly  dissimilar  on  both 
images,  then  the  probability  is  great  that  the  correlation  algorithm  will 
find  a correlation  maximum  toward  the  ends  of  the  search  segment.  Thus 
the  Am/&c  function  value  becomes  too  low  or  too  high,  resulting  in  a biased 
next  prediction.  It  has  been  observed  in  these  cases  that  succeeding  corre- 
lation maxi mums  along  a path  are  found  alternately  at  one  end  of  the  segment 
and  then  the  other,  causing  the  flu/ £pc  function  to  oscillate  rapidly  until 
good  feature  lock-on  can  be  achieved  or  until  the  prediction  mechanism 
totally  degenerates. 

For  larger  values  of  #c,  say  12  or  16  lines,*  a 7-line  search  segment 
is  unresponsive  to  rapid  terrain  changes.  In  these  cases,  the  search 
segment  length  must  be  increased.  For  the  current  matching  effort,  a 7-line 
search  segment  with  a Ax  value  of  8 lines  has  been  found  to  produce  the 
best  results.  Of  course,  this  is  very  much  dependent  upon  the  scale  of  the 
imagery  and  the  type  of  terrain  being  matched. 

3.3.3  CORRELATION  PATCH  SHAPING 

The  flu/ftc  function  is  not  only  a valid  metric  for  prediction  purposes 
but  also  for  measuring  the  amount  of  local  Image  B compression  or  expansion 
relative  to  the  same  area  on  image  A.  In  previous  block  matching  systems,  a 
correlation  search  was  performed  using  patches  of  equal  size  on  the  A and  B 
Images.  It  was  found,  however,  that  this  procedure  is  valid  only  in  flat 
terrain,  that  is,  where  tfi/flx  is  close  to  zero  and  £p / is  close  to  1.0. 

For  all  other  cases,  the  B Image  patch  must  be  compressed  or  expanded  pri- 
marily in  the  major  parallax  direction.  Referring  to  Figure  3-7  as  an  ex- 
ample, if  the  size  of  the  Image  A patch  is  chosen  as  21  x 21  cells  and  the 
current  flu/ Ax  function  is  .6,  then  the  corresponding  patch  on  the  B Image  is 


21  cells  by  13  lines.  The  Image  B patch  and  search 
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Figure  3-7  Patch  Shaping 


area  height  remains  at  21  cells  because  the  image  y parallax  is  negligible 
in  terms  of  whole  pixel  intervals. 


As  a further  sophistication  of  patch  shaping,  the  B Image  patch 
sides  (the  21-cell  dimension  of  the  figure)  may  be  slanted  to  account  for 
Image  B warp  due  to  terrain  slope  in  the  y direction. 


The  correlation  coefficient  for  these  patches  must  be  computed 
using  21  x 21  pixel  samples.  Therefore,  the  Image  B patch  area  must  be 
resampled  to  derive  21  x 21  gray-scale  values  from  the  21  x 13  pixels  that 
are  contained  in  the  patch  area.  In  the  current  Implementation  of  the 
block  matching  system,  the  entire  search  area  is  shaped  according  to  the 
£m/Ak  function  and  a synthetic  search  area  is  generated  which  contains  the 
required  number  of  samples  for  the  correlation  algorithm. 


The  need  for  patch  shaping  in  other  than  flat  terrain  is  corroborated 
In  Table  3-2.  The  table  hypothesizes  that  the  correlation  patch  size  is  the 
same  on  both  the  A and  B Images,  and  computes  the  percentage  of  unwanted 
gray-scale  samples  that  effect  the  correlation  coefficient  over  terrain 
sloping  upward  in  the  process  direction.  B Image  patch  shaping  eliminates 
these  unwanted  samples. 


average  position  of  the  blocks  In  adjacent  paths.  If  not,  the  position  of 
the  wandering  block  and  the  function  for  that  path  are  corrected 

accordingly.  Referring  to  the  figure,  the  block  on  path  j+2  has  been  de- 
tected as  wandering.  Its  position  is  corrected  to  point  P,  the  average 
poaitlon  of  paths  j+1  and  j+3. 

A special  case  exists  for  the  end  paths,  paths  j and  j+5  of  the  figure, 
where  complete  harnessing  Is  not  possible.  A problem  also  exists  when 
blocks  In  more  than  one  path  begin  to  wander.  The  probability  of  this 
occurance  on  good  imagery  is  rather  low,  but  when  it  does  occur,  the 
described  mechanism  can  correct  the  wandering  blocks  iteratively  over  a 
longer  stretch  of  imagery.  When  all  paths  become  lost,  there  is  not 
much  that  can  be  done. 


4.0  PROCESSING  RESULTS 


A digital  stereo  pair  of  Images  was  received  from  CSAETL  with  all  the 
appropriate  interior  orientation  transformations  and  exterior  orientation 
parameters.  The  digital  image  data  represents  a four  square  inch  scene  from 
the  1:48000  Phoenix-South  Mountain  model.  Digital  encoding  was  performed 
at  ETL  with  a 35  micrometer  spot  size  and  a 24  micrometer  sampling  Interval. 
A pixel,  therefore,  represents  a nominal  4 feet  on  the  ground. 


Supplied  with  the  imagery  was  also  a file  of  25,250  match  points  generated 
by  an  ETL  block  process.  These  match  points  covered  an  image  area  of  approxi- 
mately .25  inches  by  .5  inches  at  original  image  scale  and  were  used  to 
initialize  both  the  CDC  strip  process  and  block  process. 


The  modified  strip  process  and  block  process  as  described  in  the  sections 
above  were  applied  to  the  same  Image  area  and  files  of  natch  points  were 
generated.  For  the  strip  processing,  the  strip  centers  were  located  10  pixels 
apart  and  the  strip  width  was  chosen  as  20  pixels.  The  effective  recursive 
correlation  area  along  the  strips  was  chosen  as  6 scan  lines.  Even  though 
the  strip  process  provides  line-by-line  tracking  of  an  Image,  match  points 
were  output  every  8 scan  lines,  thus  providing  a match  point  lattice  whose 
Interval  is  10  pixels  in  the  y direction  by  8 lines  in  the  x direction. 


In  both  processes.  Frame  98  of  the  stereo  pair  was  chosen  as  Image  A, 
the  independent  image,  and  Frame  99  as  Image  B or  the  dependent  image.  For 
the  block  process,  a correlation  patch  size  of  21  x 21  pixels  on  Image  A 
was  chosen.  The  distance  between  blocks  or  block  paths  an  Image  A was  10 
pixels  and  /pc  or  the  block  jump-interval  was  8 lines.  A 7-line  search  seg- 
aent  was  used  on  Image  B,  including  three  scan  lines  on  cither  side  of  a - 
predicted  point. 


As  a result,  three  match  point  files  were  available  for  analysis;  the  ETL 
file,  the  block  file,  and  the  strip  file.  These  files  sere  processed  using 
photogramme trie  Intersection  to  obtain  for  each  point  at  orthometric  terrain 


elevation  in  feet  above  the  earth.  These  elevations  were  given  X,Y  coordinates 
corresponding  to  the  digital  raster  of  Image  A.  Therefore,  the  resultant 
terrain  data  is  not  truly  orthographic,  but  rather  as  unrectlfied  as  Image  A. 


These  files  of  terrain  data  were  fit  with  local,  smooth  bicubic  poly- 
nomials to  produce  the  contour  lines  shown  in  Figure  4-l(a,b,c).  The  contour 
interval  is  20  feet  and  the  background  is  the  Image  A section.  This  section 
of  image  has  been  enlarged  8 times,  the  original  scale  area  being  approximately 
•25  inches  by  .5  inches.  To  point  out  the  differences  between  approaches, 
contour  difference  images  were  generated  and  appear  in  Figure  4-l(d,e,f). 

These  difference  images  were  produced  by  subtracting  corresponding  contour 
region  images,  where  the  20  feet  contour  regions  are  alternately  colored 
black  and  white.  The  neutral  gray  color  indicates  similarity  between  contour 
regions  or  2ero  in  the  subtraction.  The  white  and  black  contour  differences 
change  polarity  on  every  contour  region.  Therefore,  visual  reference  must 
be  made  to  Figure  4-l(a,b,c)  to  determine  which  difference  in  Figure  4-l(d,e,f) 
came  from  which  image.  It  must  be  kept  in  mind  that  because  the  terrain  eleva- 
tions have  been  quantized  into  discrete  20  feet  intervals  to  produce  the 
difference  images,  the  actual  magnitudes  of  the  differences  observed  can 
range  between  0 and  20  feet.  As  a quantitative  measure,  the  RMS  magnitude  of 
the  actual  orthooetric  elevation  differences  is  in  the  neighborhood  of  8 feet 
for  all  three  cases. 


The  processing  of  the  Phoenix  model  using  the  strip  approach  and  block 
approach  was  extended  beyond  the  ETL  match  point  data  to  include  about  12,000 
match  points  across  the  entire  length  of  the  digital  sample.  The  results  of 
this  processing  in  contour  line  form  along  with  the  contour  differences  appear 
ih  Figure  4-2.  These  images  are  enlarged  4 times.  The  actual  image  area  is 
approximately  .5  by  2 inches  at  the  original  scale  of  the  photography. 
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(c)  Strip  Approach 

(f)  Block  vs.  Strip  Approach 

Figure  4-1.  Comparison  of  Hatching  Approaches  In  Terms  of  Contoured 

Terrain  Data 

a,b,c  - Contour  Images,  20  feet  Interval 

d,e,f  - Contour  Difference  Images 
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Figure  4-2.  Comparison  of  Matching  Approaches  Over  A Larger  Area 


A by-product  of  the  strip  approach  is  the  warping  of  the  B Image  to 
register  with  the  A Image  and  also  the  generation  of  a tonal  difference  image 
for  change  detection  purposes.  This  difference  image  and  the  corresponding 
unrectified  section  of  Image  B as  reference  are  shown  in  Figure  4-3. 


The  tonal  difference  image  is  a qualitative  measure  of  the  goodness  of 
match  for  the  strip  approach.  As  can  be  seen  in  the  imasc,  the  only  real  gray 
scale  differences  after  registration  are  the  reseau  marks  and  the  dissimiler 
appearance  of  the  tall  buildings  in  the  right  of  the  scene.  Referring  back 
to  Figure  4-2,  it  can  also  be  seen  that  some  of  the  buildings  were  contoured 
as  part  of  the  terrain  by  the  strip  approach.  But  despite  this  goodness  of 
registration,  the  actual  RMS  differences  in  elevation  between  the  strip 
approach  and  the  block  approach  for  the  12,000  match  points  was  again  about 
8 feet. 

As  a fcher  demonstration  of  the  performance  of  th;*:  two  approaches,  the 
raw  terrain  data,  without  fitting  or  smoothing,  that  was  produced  by  inter- 
section from  the  match  points  was  plotted  in  a three-dimensional  mode  on  a 
Calcomp  plotter.  These  plots  appear  in  Figures  4-4  and  £;-5.  Each  plotted 
line  in  these  figures  represents  one  profile  in  the  y direction.  As  was 
mentioned  before,  the  interval  between  profiles  is  8 scant  lines.  The  vertical 
dimension  in  these  plots  has  been  greatly  exaggerated  to  emphasize  small- 
scale  differences  between  the  approaches.  It  can  be  see.n  that  the  strip 
process  produces  a more  noisy  data  set. 

Regarding  processing  speed,  the  general  purpose  implementation  of  the 
block  approach  as  described  above  produces  14  match  poinits  per  central  pro- 
cessor second.  The  strip  approach  produces  27  match  poLnts  per  second.  This 
figure  for  the  strip  approach  is  based  on  a match  point  every  8 scan  lines. 

But  because  the  strip  process  performs  correlation  on  a a can  line-by-scan 
line  basis,  the  capability  exists  for  generating  a match  point  for  each  strip 
on  every  scan  line  at  no  extra  time.  Therefore,  the  process  can  produce  up  to 
216  match  points  per  second.  Both  the  block  process  ana  strip  process  are 
implemented  in  FORTRAN  running  on  the  CDC  6600  computer  system. 
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Figure  4-3.  Strip  Process  Registration  Evaluation 


5.0  CONCLUSIONS 


It  can  be  seen  from  the  contour  images  presented  above  that  both  the  atrip 
approach  and  the  block  approach  in  their  present  implementations  are  following 
the  general  terrain  trends  of  the  imagery.  It  is  conjectured  that  if  the 
terrain  data  from  the  two  approaches  were  evaluated  against  known  ground  or 
model  control  points  using  standard  photogramme trie  techniques,  the  RMS 
deviations  at  the  control  points  would  be  very  small.  This  is  due  to  the  fact 
that  control  points  are  generally  associated  with  well-defined  features,  which 
features  can  be  matched  very  well  by  an  automatic  process.  It  is  the  feature- 
less, be tween-control-point  areas  that  pose  problems  both  for  automatic 
matching  and  terrain  data  evaluation. 

The  strip  approach  is  a very  fast  production-oriented  process.  But  the 
digital  terrain  data  it  generates  is  rather  noisy.  This  seems  to  be  an 
attribute  of  tl*  line-by-line,  error  correcting  nature  of  the  process.  By- 
placing  constraints  on  the  process  smoother  terrain  data  results,  but  the 
process  becomes  less  responsive  to  terrain  fluctuations.  That  is,  the  process 
tends  to  overshoot  mountain  peaks  and  to  dig  into  valleys. 

The  block  approach,  on  the  other  hand,  generates  less  noisy  terrain  data 
due  to  its  -deterministic  predict-ahead  mechanism.  The  values  of  the  correla- 
tion coefficient  are  generally  higher  in  the  block  approach  than  in  the  strip 
approach  and  the  match  points  generated  in  feature-rich  image  areas  seem  to 
be  more  accurate  in  the  block  approach.  However,  in  featureless  and  dissimilar 
image  areas  the  block  process  has  less  system  Inertia  to  carry  it  through  the 
difficult  regions  than  the  strip  process.  In  addition,  the  block  process  is 
slower. 

• 

With  all  these  tradeoffs  in  mind,  the  block  process  is  being  considered 
further  for  implementation  in  the  digital  cartographic  benchmark. 


. 
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